/* JGrass - Free Open Source Java GIS http://www.jgrass.org 
 * (C) HydroloGIS - www.hydrologis.com 
 * 
 * This library is free software; you can redistribute it and/or modify it under
 * the terms of the GNU Library General Public License as published by the Free
 * Software Foundation; either version 2 of the License, or (at your option) any
 * later version.
 * 
 * This library is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more
 * details.
 * 
 * You should have received a copy of the GNU Library General Public License
 * along with this library; if not, write to the Free Foundation, Inc., 59
 * Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */
package org.jgrasstools.hortonmachine.modules.hydrogeomorphology.adige.hymod;

import static java.lang.Math.max;
import static java.lang.Math.min;
import static java.lang.Math.pow;
import static org.jgrasstools.gears.libs.modules.JGTConstants.isNovalue;

import java.io.IOException;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;

import org.geotools.feature.FeatureIterator;
import org.jgrasstools.gears.libs.monitor.IJGTProgressMonitor;
import org.jgrasstools.hortonmachine.modules.hydrogeomorphology.adige.Adige;
import org.jgrasstools.hortonmachine.modules.hydrogeomorphology.adige.IAdigeEngine;
import org.jgrasstools.hortonmachine.modules.hydrogeomorphology.adige.core.IDischargeContributor;
import org.jgrasstools.hortonmachine.modules.hydrogeomorphology.adige.core.IHillSlope;
import org.jgrasstools.hortonmachine.modules.hydrogeomorphology.adige.core.PfafstetterNumber;
import org.joda.time.DateTime;
import org.opengis.feature.simple.SimpleFeature;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;

/**
 * The Hymod engine for the {@link Adige} framework.
 * 
 * @author Andrea Antonello (www.hydrologis.com)
 * @author Silvia Franceschi (www.hydrologis.com)
 * @author Giuseppe Formetta
 */
public class HymodAdigeEngine implements IAdigeEngine {

	private final HymodInputs hymodInputs;
	private final List<IHillSlope> orderedHillslopes;
	// private double[] initialConditions;
	private double[] xSlow = null;
	private double[] xLoss = null;
	private double[][] xQuick = null;
	private final HashMap<Integer, double[]> outDischarge;
	private final HashMap<Integer, double[]> outVelocity;

	private final HashMap<Integer, double[]> outSubDischarge;
	private final HashMap<Integer, double[]> outDischargeInternal;
	private final HashMap<Integer, Integer> index2Basinid;
	private List<IDischargeContributor> dischargeContributorList = new ArrayList<IDischargeContributor>();
	private final boolean doPrint;
	private final boolean doLog;
	private final IJGTProgressMonitor pm;
	private final List<String> pfaffsList;
	private double[] coeffs;
	private List<Double> netnumList;
	private List<Double> distanceList;
	private HashMap<Integer, Double> netnum2Dist;
	int conta = 0;
	private int ttTimestep = 0;

	private double vel;

	public HymodAdigeEngine(HymodInputs hymodInputs,
			List<IHillSlope> orderedHillslopes,
			HashMap<Integer, Integer> index2Basinid,
			HashMap<Integer, double[]> outDischarge,
			HashMap<Integer, double[]> outSubDischarge,
			HashMap<Integer, double[]> outVelocity, List<String> pfaffsList,
			boolean doLog, boolean doPrint, IJGTProgressMonitor pm,
			double initalconditios[]) {
		this.outVelocity = outVelocity;
		this.hymodInputs = hymodInputs;
		this.orderedHillslopes = orderedHillslopes;
		this.index2Basinid = index2Basinid;
		this.outDischarge = outDischarge;
		this.outSubDischarge = outSubDischarge;
		this.pfaffsList = pfaffsList;
		this.doLog = doLog;
		this.doPrint = doPrint;
		this.pm = pm;
		outDischargeInternal = new HashMap<Integer, double[]>();

	}

	public void addDischargeContributor(
			IDischargeContributor dischargeContributor) {
		dischargeContributorList.add(dischargeContributor);
	}

	public HashMap<Integer, double[]> getDischarge() {
		return outDischarge;
	}

	public HashMap<Integer, double[]> getSubDischarge() {
		return outSubDischarge;
	}

	public HashMap<Integer, double[]> getVelocity() {
		return outVelocity;
	}

	public double[] solve(DateTime currentTimstamp, int tTimestep,
			double internalTimestepInMinutes, double[] initialConditions,
			double[] rainArray, double[] etpArray) throws IOException {
		ttTimestep = tTimestep;

		// ///////////////////////////////////////////////////////////////////////////////
		// //////////////////////////////////////////////////////////////////////////////

		netnumList = new ArrayList<Double>();
		distanceList = new ArrayList<Double>();

		/*
		 * counter for the number of station with measured value !=0.
		 */
		/*
		 * Store the station coordinates and measured data in the array.
		 */
		if (hymodInputs.pDoGeom) {
			if (hymodInputs.inSubbasinDist != null) {

				netnum2Dist = new HashMap<Integer, Double>();
				FeatureIterator<SimpleFeature> iterator = hymodInputs.inSubbasinDist
						.features();
				try {
					while (iterator.hasNext()) {
						SimpleFeature feature = iterator.next();
						int name = ((Number) feature
								.getAttribute(hymodInputs.fIdnetnum))
								.intValue();

						double z = 0;
						if (hymodInputs.fIdDistanze != null) {
							try {
								z = ((Number) feature
										.getAttribute(hymodInputs.fIdDistanze))
										.doubleValue();
							} catch (NullPointerException e) {

							}
						}
						netnum2Dist.put(name, z);
					}
				} finally {
					iterator.close();
				}
			}
		}
		// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

		if (initialConditions != null) {

			for (int i = orderedHillslopes.size() - 1; i >= 0; i--) {
				xLoss[i] = initialConditions[i];
				xSlow[i] = initialConditions[i + orderedHillslopes.size()];
				xQuick[0][i] = initialConditions[i + 2
						* orderedHillslopes.size()];
				xQuick[1][i] = initialConditions[i + 3
						* orderedHillslopes.size()];
				xQuick[2][i] = initialConditions[i + 4
						* orderedHillslopes.size()];

			}
		}

		if (initialConditions == null) {
			xSlow = new double[orderedHillslopes.size()];
			coeffs = new double[orderedHillslopes.size()];
			xQuick = new double[3][orderedHillslopes.size()];
			xLoss = new double[orderedHillslopes.size()];
			initialConditions = new double[orderedHillslopes.size() * 6];
			for (int i = orderedHillslopes.size() - 1; i >= 0; i--) {
				IHillSlope hillSlope = orderedHillslopes.get(i);
				double areaKm2 = hillSlope.getHillslopeArea() / 1E6;
				initialConditions[i + 5 * orderedHillslopes.size()] = hymodInputs.pQ0
						/ (i + 1);

				coeffs[i] = (Math.pow(10, 9)) * tTimestep * 60
						/ (areaKm2 * (Math.pow(10, 12)));
				int idd = index2Basinid.get(i);
				double RRS = 0;
				if (hymodInputs.pDoGeom) {
					RRS = hymodInputs.pVelSlow / netnum2Dist.get(idd);
				} else {
					RRS = hymodInputs.pRs;
				}

				xSlow[i] = hymodInputs.pQ0 * coeffs[i] / RRS;
				// xLoss[i]= xSlow[i]/100;
			}
		}

		for (int i = orderedHillslopes.size() - 1; i >= 0; i--) {
			IHillSlope hillSlope = orderedHillslopes.get(i);
			// /////////////////////correggereeeeeeeeeeeee////////
			double rain = rainArray[i];
			double etp = etpArray[i];

			if (rain < 0) {
				rain = 0;

			}
			if (etp < 0) {
				etp = 0;

			}

			double[] out_excess = excess(xLoss[i], rain, etp);
			double UT1 = out_excess[0];
			double UT2 = out_excess[1];
			xLoss[i] = out_excess[2];
			Integer basinId = index2Basinid.get(i);
			double UQ = hymodInputs.pAlpha * UT2 + UT1;
			double US = (1.0 - hymodInputs.pAlpha) * UT2;

			double inflow = US;
			double RRRS = 0;
			if (hymodInputs.pDoGeom) {
				RRRS = hymodInputs.pVelSlow/ netnum2Dist.get(basinId);
			} else {
				RRRS = hymodInputs.pRs;
			}
			double[] out_linres1 = linres(xSlow[i], inflow, RRRS, 1);

			xSlow[i] = out_linres1[0];
			double outflow1 = out_linres1[1];
			double QS = outflow1;
			inflow = UQ;
			double outflow2 = 0;
			double RRRQ = 0;
			if (hymodInputs.pDoGeom) {
				RRRQ = hymodInputs.pVelSlow / netnum2Dist.get(basinId);
			}
			else{
				RRRQ = hymodInputs.pRq;
			}
			for (int k = 0; k < 3; k++) {
				double[] out_linres2 = linres(xQuick[k][i], inflow, RRRQ, 1);
				xQuick[k][i] = out_linres2[0];
				outflow2 = out_linres2[1];
				inflow = outflow2;
			}

			double basinDischarge = (QS + outflow2) / (coeffs[i]);
			double basinSubDischarge = xLoss[i];

			double allContributionsDischarge = handleContributors(hillSlope,
					basinDischarge, i, initialConditions);

			/*
			 * sum together the discharge contributed by the current hillslope
			 * plus the contributions coming from upstream
			 */
			// basinDischarge = basinDischarge + allContributionsDischarge;
			basinDischarge = allContributionsDischarge;
			PfafstetterNumber pfaf = hillSlope.getPfafstetterNumber();
			if (pfaffsList.contains(pfaf.toString())) {
				outDischarge.put(basinId, new double[] { basinDischarge });
				outSubDischarge
						.put(basinId, new double[] { basinSubDischarge });
				outVelocity.put(basinId, new double[] { vel });
				// System.out.println(basinDischarge+"   sub="+basinSubDischarge+"  velocity= "+vel);

			}
			initialConditions[i] = xLoss[i];
			initialConditions[i + orderedHillslopes.size()] = xSlow[i];
			initialConditions[i + 2 * orderedHillslopes.size()] = xQuick[0][i];
			initialConditions[i + 3 * orderedHillslopes.size()] = xQuick[1][i];
			initialConditions[i + 4 * orderedHillslopes.size()] = xQuick[2][i];
			initialConditions[i + 5 * orderedHillslopes.size()] = basinDischarge;

			outDischargeInternal.put(basinId, new double[] { basinDischarge });

		}

		return initialConditions;
	}

	private double handleContributors(IHillSlope hillSlope,
			final double basinDischarge, int indice, double[] initialcond) {
		double summedContributions = 0;

		List<IHillSlope> connectedUpstreamHillSlopes = hillSlope
				.getConnectedUpstreamElements();
		hillSlope.getConnectedDownstreamElement();

		double areamonte = orderedHillslopes.get(indice).getUpstreamArea(null);
		if (connectedUpstreamHillSlopes != null) {
			for (IHillSlope tmpHillSlope : connectedUpstreamHillSlopes) {
				PfafstetterNumber pNum = tmpHillSlope.getPfafstetterNumber();
				int hillslopeId = tmpHillSlope.getHillslopeId();

				double upstreamDischarge = outDischargeInternal
						.get(hillslopeId)[0];

				/*
				 * handle the contributors
				 */
				for (IDischargeContributor dContributor : dischargeContributorList) {
					Double contributedDischarge = dContributor
							.getDischarge(pNum.toString());
					if (!isNovalue(contributedDischarge)) {
						if (doLog && doPrint) {
							pm.message("----> For hillslope "
									+ hillSlope.getPfafstetterNumber()
									+ " using hydrometer/dams data in pfafstetter: "
									+ pNum.toString() + "(meaning added "
									+ contributedDischarge + " instead of "
									+ upstreamDischarge + ")");
						}

						/*
						 * here the contributor will give its contribution,
						 * which depends on the type of contributor. For example
						 * a Hydrometer will completely substitute the
						 * calculated discharge of the current hillslope
						 * (tmpHillSlope) with the measure supplied by the
						 * Hydrometer.
						 */

						// funziona
						// if (contributedDischarge != -9999) {
						// upstreamDischarge = dContributor
						// .mergeWithDischarge(contributedDischarge,
						// upstreamDischarge);
						// }

						// modificato
						if (!isNovalue(contributedDischarge)) {

							upstreamDischarge = dContributor
									.mergeWithDischarge(contributedDischarge,
											upstreamDischarge);
						}

					}
				}

				summedContributions = summedContributions + upstreamDischarge;

			}
		}
		double output = 0;
		boolean pdoRouting = hymodInputs.pDorouting;
		if (pdoRouting == true) {
			output = doRouting(summedContributions, basinDischarge, hillSlope,
					initialcond[indice + 5 * orderedHillslopes.size()],
					areamonte);
		} else {
			output = basinDischarge + summedContributions;

		}

		return output;
	}

	private double doRouting(double UPdischarge, double HYMODbasinDischarge,
			IHillSlope hillslope, double condini, double Areaamonte) {

		Areaamonte = hillslope.getHillslopeArea();
		double a = 0;
		double b = (ttTimestep * 60);

		int steps = 50;
		double h = (b - a) / steps;
		double k1 = 0, k2 = 0, k3 = 0, k4 = 0;
		double y;
		int i;

		y = condini;

		double lambda1 = 0, lambda2 = 0, Qr = 0, vr = 0;
		if (hymodInputs.LAMBDA1 != null) {
			lambda1 = hymodInputs.LAMBDA1;
		} else {
			lambda1 = 0.2;

		}
		if (hymodInputs.LAMBDA2 != null) {
			lambda2 = hymodInputs.LAMBDA2;
		} else {
			lambda2 = -0.1;

		}
		if (hymodInputs.QR != null) {
			Qr = hymodInputs.QR;
		} else {
			Qr = 100;

		}
		if (hymodInputs.vR != null) {
			vr = hymodInputs.vR;
		} else {
			vr = 1.5;

		}

		double Areatotale = orderedHillslopes.get(0).getUpstreamArea(null);
		Areatotale = 1;
		double l = hillslope.getLinkLength();
		double deriv = 0;

		for (i = 0; i < steps; i++) {

			double TTT = T(hillslope, Areaamonte, lambda1, lambda2, vr, Qr,
					Areatotale, l, y);

			deriv = -TTT * y + TTT * HYMODbasinDischarge + TTT * UPdischarge;

			k1 = h * deriv;

			deriv = -TTT * (y + k1 / 2.0) + TTT * HYMODbasinDischarge + TTT
					* UPdischarge;

			k2 = h * deriv;

			deriv = -TTT * (y + k2 / 2.0) + TTT * HYMODbasinDischarge + TTT
					* UPdischarge;

			k3 = h * deriv;

			deriv = -TTT * (y + k3) + TTT * HYMODbasinDischarge + TTT
					* UPdischarge;

			k4 = h * deriv;

			y += k1 / 6.0 + k2 / 3.0 + k3 / 3.0 + k4 / 6.0;

		}

		vel = vr * (Math.pow(y / Qr, lambda1))
				* (Math.pow(Areaamonte / Areatotale, lambda2));

		return y;

	}

	private double T(IHillSlope hillslope, double areaamonte, double lambda1,
			double lambda2, double vr, double Qr, double Ar, double l,
			double portata) {
		double re = 0;

		double vvv = vr * Math.pow(portata, lambda1)
				* Math.pow(hillslope.getHillslopeArea(), lambda2);
		double kkk = vvv / (l * (1 - lambda1));
		if (kkk == 0) {
			re = 0.00001;
		} else {
			re = kkk;
		}
		return re;

	}

	private double[] excess(double x_losss, double Pval, double PETval) {
		double[] o_exces = new double[3];
		double pB = hymodInputs.pB;
		double pCmax = hymodInputs.pCmax;

		double xn_prev = x_losss;
		double coeff1 = ((1.0 - ((pB + 1.0) * (xn_prev) / pCmax)));
		double exp = 1.0 / (pB + 1.0);
		double ct_prev = pCmax * (1.0 - pow(coeff1, exp));
		double UT1 = max((Pval - pCmax + ct_prev), 0.0);
		Pval = Pval - UT1;
		double dummy = min(((ct_prev + Pval) / pCmax), 1.0);
		double coeff2 = (1.0 - dummy);
		double exp2 = (pB + 1.0);
		double xn = (pCmax / (pB + 1.0)) * (1.0 - (pow(coeff2, exp2)));
		double UT2 = max(Pval - (xn - xn_prev), 0);
		double evap = min(xn, PETval);

		xn = xn - evap;

		o_exces[0] = UT1;
		o_exces[1] = UT2;
		o_exces[2] = xn;
		if (xn != xn || UT1 != UT1 || UT2 != UT2) {
			System.out.println("ATTENTION: NaN for some system variable");

		}

		return o_exces;

	}

	private double[] linres(double x_sloww, double infloww, double RR, double dt) {

		double[] o_linres = new double[2];

		double x_sloww_prev = x_sloww;
		double x_slow_new = (infloww * (1.0 - RR)) + x_sloww_prev * (1.0 - RR);
		double outfloww = x_slow_new * RR / (1 - RR);

		o_linres[0] = x_slow_new;
		o_linres[1] = outfloww;
		return o_linres;

	}

}
